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ABSTRACT 

Since it has become possible to discovery planets orbiting nearby solar-type 
stars through very precise Doppler-shift measurements, the role of methods 
used to analyze such observations has grown significantly. The widely employed 
model-dependent approach based on the least-squares fit of the Keplerian motion 
to the radial-velocity variations can be, as we show, unsatisfactory. Thus, in this 
paper, we propose a new method that may be easily and successfully applied 
to the Doppler-shift measurements. This method allows us to analyze the data 
without assuming any specific model and yet to extract all significant features 
of the observations. This very simple idea, based on the subsequent subtraction 
of all harmonic components from the data, can be easily implemented. We 
show that our method can be used to analyze real 16 Cygni B Doppler-shift 
observations with a surprising but correct result which is substantially different 
from that based on the least-squares fit of a Keplerian orbit. Namely, using 
frequency analysis we show that with the current accuracy of this star's 
observations it is not possible to determine the value of the orbital eccentricity 
which is claimed to be as high as 0.6. 

Subject headings: Stars: Individual (16 Cyg B) - Planetary Systems, methods: 
numerical 
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1. Introduction 

Recent improvements in the long-term precision of Doppler-shift measurements (Butler 
et al. 1996) resulted in several spectacular detections of planetary companions to solar-type 
stars (for review see the paper of Marcy & Butler 1997). As such discoveries supply 
indirect evidence of the existence of extra-solar planets, other explanations of observed 
radial- velocity variations appeared, e.g., stellar pulsations (Gray 1997). The most recent 
results, however, show that only a planetary hypothesis is acceptable (Marcy 1998; Gray 
1998). 

The usual procedure showing that there exists a planet around s star consits of a 
direct least-square fit of Keplerian model to the observations. It always gives certain 
values of orbital parameters and their formal errors. In the case of 'good' data this is the 
best and the quickest way to obtain relable results. However, in the case of spare data 
with big errors one has to prove that the least-squares method can be used and that the 
obtained parameter values and their errors are good estimates of the real values. This is a 
difficult and time consuming task. Without doubt, we have this situation with the Doppler 
observation of extra-solar planets. We present an analysis of this problem and we show that 
the eccentricity of the fitted orbit is a very sensitive parameter and, in some cases, its value 
and error given by the least-squares method are not correct. 

The aim of this paper is to show how the mentioned problem can be solved in 
practice. Namely, we propose a method that can be very useful for analyzing radial-velocity 
variations. It is based on a simple idea involving the subsequent subtraction of periodic 
components from the data. This approach allows us to analyze the observations without 
assuming any specific model describing the system behavior (like Keplerian motion or 
stellar pulsations). After the determination of all significant components of the data, it 
remains to be decided which process is responsible for what we observe and whether it is 
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possible to choose only one. 

The plan of this paper is as follows. In Section 2 we analyse observations of 16 Cygni 
B and we explain why the standard least-squere fit does not give reliable estimates of 
parameters and their errors. In Section 3 we analyze analytically the Keplerian motion of 
the system 'a star with one planet' in order to learn how its motion modulates the observed 
star radial-velocities. We investigate mainly the spectral properties of the motion which are 
essential for our method. In Section 4 we develop a simple numerical technique which can 
be used to extract all the information we need to compare with the results from Section 
3. In Section 5, we perform a numerical test of the method using simulated radial-velocity 
variations with the orbital parameters of 70 Vir (Marcy & Butler 1996). In section 6, we 
discuss the application of the method to finding the eccentricity of 16 Cygni B (Cochran 
et al. 1997). 

2. Least-squares analysis of radial velocities of 16 Cygni B 

In order to determine the quality and reliability of the least-squares fit of a Keplerian 
orbit to the radial velocity measurements of 16 Cygni B performed by Cochran et al. 
(1997) we analyzed the topology of the x 2 minimum on the (a, e) plane, where a is the 
semi-major axis and e is the eccentricity. To this end, we used the real observations 
published in Cochran et al. (1997) and the Levenberg-Marquard method to solve the 
nonlinear least-squares problem (in our case a fit of a Keplerian orbit) which cannot be 
linearized (as this is the case for 16 Cygni B observations, as we show below). We took the 
necessary FORTRAN code from MINPACK library (More et al. 1980). As the semi-major 
axis and the eccentricity are the most crucial parameters of the model, we compute the 
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behavior of 

2 \ modeled observed J 

~i 

on the (e, a) plane. For a given point (a, e) other parameters (P, T p , oj) are always chosen 
to correspond to the global minimum of x 2 , i-e., for fixed values of parameters (a, e) we 
make a series of fits with initial values of the elements (P, T p , u) covering their whole range 
and we take the smallest value of x 2 ■ 111 this way we get the behavior of x 2 on the (e, a) 
plane, see Figure 1. As one can see the confidence levels la, 2a, 3a of the parameters (e, a) 
bound a large region of the parameters' plane. This shows that e is very weakly determined 
with the data available for 16 Cygni B. At this point we make a well known but important 
remark. Namely, as one can see in Figure 1, the confidence lines do not determine ellipses 
as it should be if we could use la, 2a, 3a levels to determine the errors of a and e (Press 
et al. 1992). It also means that the errors of the parameters obtained from the linearized 
(or not) least-squares fit cannot be treated as the correct estimates of the parameters' 
accuracy. Thus, for example, the value 0.082 from the paper of Cochran et al. (1997) 
has little to do with the accuracy of the determined eccentricity. In fact, using nonlinear 
Levenberg-Marquard least-squares method we find a better solution for the same data of 16 
Cygni B (this solution is marked with a filled square in Figure 1; the filled star indicates the 
orbital solution found by Cochran et al. (1997)). One can compare the x 2 behavior with 
a very similar picture obtained by means of the bootstrap method (Press et al. 1992), see 
Figure 2. Clearly the distribution of a and e is not normal. Moreover, from this distribution 
we learn that any value of the eccentricity within the interval approximately (0.25, 1.0) is 
probable at the 3a level. 

To sum up, one can always use the least-squares method to obtain the orbital 
parameters. This method always gives a solution but sometimes the obtained parameters 
together with its errors have little to do with the real ones, whatever they are. However, 
the process of checking if the obtained result is correct is time consuming, and thus it seems 
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resonable to find a simple method which can give the correct answer quickly. We address 
this problem in Sections 3 and 4. 



3. Theoretical Background of Frequency Analysis 

Let us assume that we have a planetary system consisting of a star and one planet. 
Choosing the reference frame placed in the center of mass of this system (barycentric 
system) with the Z-th axis directed from the observer, we can calculate the Z-th coordinate 
of the star from the following equation 

ut) = ~m (i) 

where m* is the mass of the star, m is the mass of the planet and Z is its Z-th. coordinate. 
Next, because of the Keplerian motion of the system, we have 

Z(t) = a sin i (sin a; (cosE — e) + cosc^Vl — e 2 sin E^j , (2) 

where a,i,u,e are Keplerian elements of the planet (semi-major axis, inclination, longitude 
of periastron and eccentricity) and E is its eccentric anomaly that can be calculated from 
the Kepler equation 

2tt 

E -esinE = M, M = n(t-T p ), n = — , (3) 

where M is the mean anomaly, P is the orbital period of the planet and T p is the time 
of periastron. It was shown by Konacki & Maciejewski (1996) that function (1) can be 
expanded in the following series 

Z*{t) = -— (z + *ff Z k e^-^) , (4) 

where 



Zk = 2l ^ Jk -^ ~ -WfeO] + iS[J k ^(ke) + J k+1 (ke)}}e- iknT * 
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and 



3 

Z = — easinisinu;, C = a sin i sin u;, S = —a\J 1 — e 2 sin i cos u>. 
2 



In the above J n (z) is a Bessel function of argument z. 

This expansion has a very important property — terms corresponding to successive 
harmonics have decreasing amplitudes. In this way the term with the frequency n has a 
larger amplitude than that with frequency 2n, which has larger amplitude than that with 
frequency 3n, etc. Generally, it is possible to prove the following inequality (Konacki & 
Maciejewski 1996) 



MZ) = ^ = ^^ k (e,u)<l, (5) 



for each e G (0, 1) and uj G [0, 2tt); in (5) we denote 



pn^u[J' k+l ((k+ l)e)] 2 + [J k+1 ((k + l)e) 



12 



(6) 



\ (3H&n 2 u[J' k {ke)] 2 + [J k (ke)} 

and (3 2 = e 2 /(l — e 2 ). When e = only the first term with frequency n has a non-zero 
amplitude. 

These considerations allow us to state the following: Keplerian motion of the planet 
orbiting a star leads to specific changes in the Z-ih coordinate of a star. These changes have 
a very characteristic spectrum in which the term with frequency equal to the mean motion 
of the planet n is dominant and amplitudes of subsequent harmonics (with frequencies 
2n, 3n, . . .) decrease strictly monotonically. 

The expansion of the Keplerian motion we show above can be easily applied to calculate 
the radial velocity of the star resulting from its motion in the system. Having the expansion 
of the Z-ih coordinate of the star, we derive its radial-velocity by differentiating (4) with 
respect to time 

vrit) = = n— *ff W fe e^-H (7) 

x h— — oo 
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As it can be easily noticed, the above expansion has the same feature as the expansion of 
the Z-th coordinate. In fact, the ratio Ak(v r ) of two successive harmonics is always less 
than 1, since we have (Konacki & Maciejewski 1996) 

am = |Zfc ^n| 1)n| = nr Mz) = * fe(e ' u) < L (8) 

Thus, all useful properties of the radial velocity expansion can be derived from that of the 
Z-th coordinate. 

In summary, if we observe radial-velocity variations that are of planetary origin then in 
their spectra we will detect the basic frequency corresponding to the planet orbital period 
and its harmonics with monotonically decreasing amplitudes. The ratios of successive 
harmonics depend on the value of the eccentricity and the longitude of periastron. It means 
that in the case of a circular orbit we will be able to notice only one periodic term (which is 
obvious) and with increasing eccentricity the number of detectable harmonics will increase. 
Finally, we should mention that, because of the finite accuracy of our observations, we can 
only detect a few higher harmonics. Thus, from the observational point of view, expansion 
(7) is always finite and includes the basic term (corresponding to the planet orbital period) 
and a few of its harmonics. 

Let us note that we can examine characteristic features of the spectra of other processes 
(e.g. stellar pulsations) and compare them with the expansion of Keplerian induced 
radial-velocity variations. Such information might be crucial for the proper interpretation 
of observations. 

4. Numerical Realization of Frequency Analysis 

Let V° denote a set of radial velocities of a star obtained from the Doppler- 
shift measurements. Using the least-squares method, we fit the function 
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Fi(t) = Ai sin(27r/i(/: — 1 )) + B 1 sin(27r/ 1 (t — 1 )) to the data. As the first approximation of 
fi, we take the frequency corresponding to the maximum of the Lomb-Scargle periodogram 
of V° (Lomb 1976; Scargle 1982; Press et al. 1992). 

At the k-ih stage of this algorithm, we have a set of the residuals V k obtained by fitting 
function F k to the original observations, where 

F k (t) = F fc _i(*) + A k sin(27r/ fe (t - t )) + B k cos(2tt f k (t -t )). (9) 

This means that at the k-th stage of the algorithm we have the following model function 
F k (t): 

k 

F k (t) = Aj sm(2ir fj(t -t )) + B j cos(2tt fj(t - t )). (10) 
i=i 

Then, using the periodogram of V fc , we approximate fk+i, and we fit function F k+ i to the 
original set of radial-velocities. These steps are repeated until the desired number of terms 
is obtained or until the final residuals are smaller than the assumed limit. We call the above 
algorithm unconstrained frequency analysis, as we assume that all frequencies we find are 
independent. However, we can modify our method to the constrained form. To this end, we 
assume that all terms have frequencies that are natural combinations of the chosen basic 
frequency /& 

fi=jf b , ./ 1 /■•• (11) 

Thus the model we fit to the observations has the form 

k 

v r (t) = Aj sm(2njf b (t - t„)) + B 3 cos(2nj f b (t - t )). (12) 
i=i 

In this way parameters of the model are f b , Aj, Bj and not fj,Aj,Bj] j = 1, . . . , k. We call 
this version of the algorithm the constrained frequency analysis. 

Numerical tests of this method can be found in (Konacki & Maciejewski 1996) and its 
successful applications to the PSR B1257+12 pulsar timing observations in (Maciejewski & 
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Konacki 1997; Konacki et al. 1998). In the next section we show how this method works on 
a simulated set of radial velocity measurements. 

5. An Example - 70 Virginis 

For the purposes of this example, we have chosen the orbital parameters of 70 Vir 
(Marcy & Butler 1996). They are presented in Table 1 where we also present several 
periods and corresponding amplitudes calculated from expansion (7). For our tests, we 
computed data that are very similar in nature to the real observations. Thus we have 39 
unevenly sampled radial-velocities. We also added 10 ms" 1 Gaussian noise resembling the 
real observational error. 

In Figure 3 we present results from the unconstrained frequency analysis for the data 
described above. What can we learn about the planetary system from these results? First, 
we notice that the eccentricity of the orbit of 70 Vir is large since we are able to detect 
three periodic terms in the data. In fact, from the ratio of the amplitudes of the basic term 
to its first harmonic, we can precisely calculate the eccentricity (see Konacki & Maciejewski 
1996; Konacki et al. 1998). 

Let us assume now that we do not know that these radial-velocity variations are of 
planetary origin. From the unconstrained frequency analysis for 70 Vir we find out that 
there are three periodic terms detectable, and, up to the accuracy of the method, they are 
all natural combinations of a certain basic frequency. This basic frequency corresponds 
to the term with the largest amplitude. It means that in fact, we see effects from the 
basic frequency and its two subsequent harmonics. Moreover, the ratio of their amplitudes 
resembles a process of planetary (Keplerian) origin. It is extremely unlikely that there are 
other reasonable processes that can produce such spectra and in this way mimic a planet. 
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Thus, taking into account all the facts, the most natural answer would be that we observe 
radial- velocity variations induced by a planet in an eccentric orbit. 

In summary, frequency analysis can supply us with strong constraints on the possible 
models of observed radial- velocity variations and verify the orbital solution obtained from a 
least-squares fit. We should also notice that we can apply the constrained frequency analysis 
as well; it usually results in a better fit and can be used to obtain better approximations 
of values of the amplitudes and basic frequency. Finally, one can compare this result with 
the behavior of x 2 i n Figure 4. The comparison shows that in many cases the frequency 
analysis will not be superior to the least-squares fit but will give an idependent insight 
into the data. Therefore, it is appropriate to use both methods while analyzing radial 
velocity measuerements, as the frequency analysis is numerically not demanding and has 
the capablility of confirming the least-squares' findings. 

We should also mention that it is possibible to interpret the observations in a different 
way which might be as well justified as the original assumption. Namely, it is possible 
to show that in some cases a planetary system consisting of two (or more) planets can 
be misinterpreted as one planet in a highly eccentric orbit and vice versa. This is easily 
understandable since from the frequency analysis we know that an eccentric orbit just 
produces a finite set of harmonics in the real data which subsequently might be interpreted 
in different manners. As there are some difficulties with the explanation of the existence of 
planets with highly eccentric orbits (see the paper of Marcy & Butler (1998) and reference 
therein for a discussion of this topic) and there is no obvious reason why the planets 
discovered around the Solar-type stars should be the only ones in their systems the above 
idea might be worth checking. 



6. Frequency Analysis of 16 Cygni B observations 



We applied the frequency analysis to real observations of 16 Cygni B. According to 
Cochran et al. (1997) the planet is believed to move in an highly eccentric orbit (e=0.634). 
We calculated theoretical amplitudes of the first harmonics of expansion (7), presented in 
Table 2. For calculation we took the orbital parameters obtained by means of the least 
squares method by Cochran et al. (1997). Comparing the values of amplitudes from 
Table 2 and the mean error of observations we concluded that we could detect at least 
the main frequency n and its harmonic 2n. We performed the frequency analysis of the 
real observations and we were only able to detect the basic frequency corresponding to 
the orbital period (see Figure 5). After the subtraction of this term, there were no other 
dominant frequencies present in the data. Thus, our result obtained with the help of 
frequency analysis was inconsistent with result of Cochran et al. (1997) who reported the 
eccentricity 0.634 with uncertainty 0.082. This is, however, consistent with the results from 
Section 2. Stricly speaking, the frequency analysis reveals that the current quality of the 
data available for 16 Cygni B does not allow us to determine the eccentricity of the orbit 
due to the absence of any harmonics and confirms the findings based on the x 2 topology 
and the bootstrap method for these observations. 

The final conclusion is that the real error of the eccentricity is much greater than that 
reported by the least-squares method, making the determination of the eccentricity of 16 
Cygni B almost impossible. This fact can be easily derived from the frequency analysis. The 
absence of any harmonics just indicates that the eccentricity must remain undetermined. 
With the frequency analysis one gets this correct result without the time consuming x 2 
map or the bootstrap method. 



7. Conclusions 



In this paper we have shown that results obtained from least-squares fits of Keplerian 
orbits to real Doppler-shift measurements may lead to incorrect interpretations. Specifically, 
they may give unrealistic or even entirely false values of parameters and their uncertainities. 
In order to solve these problems we have proposed a new method, frequency analysis, which 
efficiently provides an independent test of the reliability of determined orbital parameters. 
This method may deliver a substantial revision of the current values of planets' high 
eccentricities that are essential for our understanding of the formation and evolution of 
planetary systems. It might even lead to hints that some of the observed high eccentric 
planets are in fact planetary systems consisting of more than one planet or at least provide 
an independent point of view on the same data. These facts, together with the ease of 
applicability of frequency analysis, make our method worth trying on future observations if 
not for the data already gathered. 



This work was supported by KBN grant 2P03D.023.10. 
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Fig. 1. — x 2 ma P 011 the plane (e, a) of the real data of 16 Cygni B. The filled square indicates 
the global minimum that is different than the orbital solution found by Cochran et al. (1997) 
denoted by the filled star. The significance levels "la", "2a" and "3a" are presented. 

Fig. 2. — Bootstrap estimation of the distribution of (e, a) based on the sample consisting of 
10 4 synthetic sets of data (left). Distribution of e (right, top) and a (right, bottom) calculated 
from the same sample. The distributions for both parameters are clearly not normal. 

Fig. 3. — Unconstrained frequency analysis of the fake reflex velocity of 70 Vir. The 
simulated data consist of 39 points from 1988 to 1996.25, 10 ms~l Gaussian noise was added. 
a, b, c, d: Subsequent steps of the frequency analysis. Found periods and their amplitudes are 
presented. Numbers in parentheses are 3a uncertainties of the parameters. There are three 
detectable periodicities in the data. 

Fig. 4. — x 2 ma P 011 the plane (e, a) of the fake data of 70 Virginis. The filled square 
indicates the global minimum found and the star indicates the parameters assumed while 
calculating the fake data. The significance level 3a is presented. 

Fig. 5. — Unconstrained frequency analysis of the real reflex velocity of 16 Cygni B. The 
data taken from the paper of Cochran et al. (1997) consist of 70 points from about 1988 to 
1997. a, b: Two subsequent steps of the frequency analysis. As one can see, there is only 
one periodic term detectable (a) corresponding to the planetary period of about 800 days. 
After subtracting this periodicity there are no significant peaks in the data present. Solid 
and dash-dotted line correspond to 90 an 50 percent significance levels, respectively. 
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Table 1. Orbital and Theoretical Frequency Analysis Parameters of 70 Vir 



70 Vir 



P (days) 


116.67 


k = 


1 


PI 


= 116 d .67 


A P1 = 


261.36 ms" 1 


T p (JD) 


2448990.403 


k = 


2 


P2 


= 58 d .335 


A P2 = 


101.08 ms- 1 


e 


0.40 


k = 


3 


P3 


= 38 d 89 


A P3 = 


43.83 ms" 1 


uj (deg) 


2.1 


k = 


4 


PA 


= 29 d .1675 


A P4 = 


19.99 ms" 1 


ai sin % (AU) 


0.00312 


k = 


5 


P5 


= 23 d 334 


A P5 = 


9.39 ms" 1 



Table 2. Orbital and Theoretical Frequency Analysis Parameters of 16 Cyg B 



16 Cyg B 



P (days) 


800.8 


k = 


1 


PI 


= 800 d 8 


A P1 


= 28.89 ms" 1 


T p (JD) 


2448935.3 


k = 


2 


P2 


= 400 d .4 


A P2 


= 16.14 ms" 1 


e 


0.634 


k = 


3 


P3 


= 200 d .2 


A P3 


= 10.23 ms- 1 


uj (deg) 


83.2 














K (ms" 1 ) 


43.9 
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Fig. 4.- 
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